Phase locked solutions and their stability in the presence of propagation delays 
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We investigate phase-locked solutions of a continuum field of nonlocally coupled identical phase 
oscillators with distance-dependent propagation delays. Equilibrium relations for both synchronous 
and traveling wave solutions in the parameter space characterizing the non-locality and time delay 
are delineated. For the synchronous states a comprehensive stability diagram is presented that 
provides a heuristic synchronization condition as well as an analytic relation for the marginal stability 
curve. The relation yields simple stability expressions in the limiting cases of local and global 
coupling of the phase oscillators. 
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INTRODUCTION 



An ensemble of coupled phase oscillators provides a simple yet powerful paradigm for studying the collective 
behaviour of many complex real life systems and has been extensively studied in this context [1-3]. Under appropriate 
conditions, phase oscillators can exhibit stable in-phase synchronous or other type of phase- locked solutions [4]. In a 
class of discrete networks of oscillators it has been shown that in-phase synchronous behavior can be achieved even 
in the presence of a fixed signal transmission delay [5] . More intricate dynamics have been discovered in nonlocally 
coupled continuous networks, whereby phase-locked and incoherent activity can simultaneously exist at different spa- 
tial locations [6], giving rise to a spatio-temporal pattern that has been termed as a chimera state [7]. In this paper, 
we discuss the existence and stability of phase-locked solutions in a continuum of nonlocally coupled identical phase 
oscillators with distance-dependent delays. Such delays are a natural consequence of the finite speed of information 
propagation in space. In similar contexts, distance-dependent delays have been considered in continuum models of 
neural activity [8, 9]. In particular, the effects of distributed delays can be very different from fixed delays [9-12], and 
the analysis is generally more difficult. 



MODEL SYSTEM AND ITS SYNCHRONOUS STATES 



We consider a continuum of identical phase oscillators, arranged on a circular ring C and labeled by x G 
with periodic boundary conditions, whose dynamics is governed by 
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where 4'{x, t) € [0, 27r) is the phase of the oscillator at location x and time t, whose intrinsic oscillation frequency 
is w > 0, if is the coupling strength and G : [—L, L] — > M is an even function describing the coupling kernel. The 
quantity v denotes the signal propagation speed which gives rise to a time delay of \z\/v for distance |z| from the 
locations x. As the oscillators are located on a ring with circumference 2L, the distance between any two oscillators 
is given by the shorter of the Euclidean distance between them along the ring. In this configuration, the maximum 
distance between the coupled oscillators is L and thus the maximum time delay would he Tm = L/v. If we denote the 
location of any two oscillators as x and x' then z = x — x' . a is a parameter which makes the coupling phase-shifted 
and has been crucial for observing chimera states. 

We choose the coupling kernel G{z) to have an exponentially decaying nature and its normalized form is taken to 

be 
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where ct > is the inverse of the interaction scale length and is a measure of the nonlocality of the coupling. The 
exponential form of G{z) and the sin coupling function are the natural consequences of the reduction of a more general 
reaction-diffusion system to a phase model under nonlocal and weak coupling limits [6]. We further make time and 
space dimensionless in Eqs.(l) and (2) by the transformations t — )■ Kt, uj ^ ur/K, k — crL, z — > z/L, — > KTm and 
X — > x/i and obtain 



d 

— (j)(x, t) = u} — / Giz) sin [Mx, t) ~ (j) (x — z,t — \z\t„i) + a] dz. 



(3) 



G(z) = -e-'^l^l (4) 

^ ' 2(1.0-6-'=) ^ ' 

We look for phase-locked solutions of Eq. (3) that have the form: 

(jyn.kix^t) = Qt + nkx + (t>o. (5) 

These solutions are phase-locked, in the sense that the difference of the phases at two fixed locations in space does 
not change in time [4]. They can describe spatially uniform equilibria (fi = fc = 0), spatial patterns (fi = 0, fc 7^ 0), 
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synchronous oscillations (fi 7^ 0, A; = 0), or traveling waves (Jl ^ 0, fc ^ 0). Notice that k needs to be an integer 
because of periodic boundary conditions, and the value of 00 can be taken to be zero by a translation. 

It is useful to first look at the undelayed case (w = 00). Substitution of (5) into (3) gives the relation between the 
temporal frequency Q and the wave number k as 

Q = OJ — J G{z) sin (Trfcz + a) dz 

= w — G(fc)sina (6) 

where G{k) = J'^ G{z) cos(7rfcz) dz denotes the discrete cosine transform of G. Note that the role of uj is simply to 
shift the value of $7, so it can be taken to be zero without loss of generality (we will see later that this is no longer 
true in the presence of delays). When a = 0, f2 is identical to the individual oscillator frequency w. More generally, 
for a € M, (6) yields a unique = Vl{k) for any given k G Z. Hence there exists a countably infinite set of solutions 
0n(fc),fei fc S Z, of (3). The condition for spatial patterns with wave number k is 

UJ = G{k) sin a, (7) 

and can only be satisfied if | sina| is not too small. Generically, however, the value of from (6) is nonzero, so the set 
{4>Q.(k),k '■ k G Z} includes a unique synchronous solution {k = 0) and the rest correspond to traveling waves (/c ^ 0), 
with wave speed equal to fl{k)/k. It turns out that in general only a few of the solutions (j)n^k can be stable. The 
linear stability is determined by the variational equation 

d 

—u{x,t) = — / G(z) cos (vrfcz + a) [u{x,t) — u{x — z,t)]dz 
where u{x,t) = (f>{x,t) — (l)Qj^{x,t). With the ansatz u{x,t) = ^^t^nrnx^ A G C, n G Z, we have 

G(z)cos(7rfcz + a)(l-e-*""^)(iz. (8) 
After some simplification and using the fact that G is an even function, one obtains 

R \ f ^ G(fc + n) + G(fc-n)\ 
ReA = —(cos a) I G(k) 1 . 

Note that the pair (A,n) ~ (0,0) is always a solution, corresponding to the rotational symmetry of the solutions (5). 
Hence, 0o_fe is linearly asymptotically stable if and only if the right hand side of (8) is negative for all nonzero integers 
n. For instance, if G has a global maximum (respectively, minimum) at fc, then the corresponding solution (t>Q,k is 
stable provided cos a > (resp., < 0). Furthermore, if G(fc) > 4|G(n)| for all n 7^ fc, then (f)i},k is the only stable phase- 
locked solution (up to a phase shift) for cos a > 0. In particular, if the coupling kernel is constant (that is, G(fc) > 
iff fc = 0) or has a predominant constant component (G(0) > 4|G(n)|, V n 7^ 0), then the synchronous solution (f>n{o).o 
is the only stable phase-locked solution (up to a phase shift) for cos a > 0, which explains the widespread use of 
globally coupled oscillators in synchronization studies. 

We now show that the inclusion of propagation delays offers a much richer solution structure. We again seek 
solutions 4>Q^k of the form (5) when the propagation velocity v is finite. The relation (6) between the temporal 
frequency and the wave number now takes the form 



n = OJ — J G{z) sin [riTm|z| -|- T^kz -\- a] dz, (9) 

which is an implicit equation in fl. Note that the right hand side is a bounded and continuous function of $7; therefore, 
(9) has a solution f2 for each A: G Z. However, in contrast to the undelayed case, it is now possible to have several 
solutions n for a given wave number fc. (Note, however, that the condition for spatial patterns (7) remains identical 
to the undelayed case.) Furthermore, the linear stability of the phase-locked solutions 4>n^k is determined through a 
dispersion relation 



^1 

A = - 



G(2) cos (rjT„,.|z| +7rfcz + a) (l-e-^l^l^'-e"""^) dz. (10) 
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FIG. 1. Frequency f7 of the phase-locked solutions given by of Eq. (9) is plotted as a function of the mean time delay r for 
K — 2.0 and a = 0.0. The panel (a) on the left is for the travelling wave number fc = 1 whereas the panel (b) on right is for 
synchronous solutions with k = 0. The different curves correspond to different values of the intrinsic oscillator frequency cj. 
The solid portions denote stable states and the dotted ones (in red) unstable states in panel (b). 



As before, linear stability requires Re(A) < for all solutions A of (10) for all nonzero integers n. The main difference 
with the undelayed case (8) is that (10) is an implicit equation in A, so its solution is not straightforward. Furthermore, 
even for a fixed value of n g Z, (10) generally has an infinite number of solutions for A. For simplicity we take a = 
and carry out the integration in equation (9) for fc = to get an analytical equation for the synchronous solutions 
as: 

Q = uj — J G{z) sin (r2r,„|z|) dz 

flTm — flTmE^'^ COs{flTm) — He~'^ sin(f2Tm) 
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which is an implicit equation in 51. Being a transcendental equation its solution can in principle be multi-valued in 
for a given set of parameters a;,T„j and k and can lead to higher branches of collective frequencies as pointed out 
by Schuster and Wagner [13] for a system of two coupled oscillators. A similar (much lengthier) analytical expression 
can also be obtained for finite values of k which we do not write here. We further define a mean delay parameter by 

1 

G{z)T^\z\dz (12) 

-1 

which weights the individual delays with the corresponding connection weights. With the exponential connectivity 
given by Eq. (4), this translates into 

T = T„— — . (13) 

K[e^ — 1) 

This gives values for the limiting cases: f — for local {k — >■ oo) and f Tm/2 for global (k — > 0) coupling. 

Fig. 1(a) and Fig 1(b) show plots of the numerical solutions of the equilibrium relations for as a function of f for 
traveling wave with k = land for synchronous {k = 0) states respectively. In both cases the value of k is taken to be 
2.0 and the plots are obtained for several values of oj. In these figures as the curves for w = 0.8 show, it is possible to 
have multiple solutions f2 for a given value of f . The stability of these higher states will be discussed in later sections 
of the paper. One also notes that the lowest branch shows frequency suppression as a function of the mean time delay 
f. 
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STABILITY OF THE SYNCHRONOUS SOLUTIONS 



We now examine the stability of the synchronous solutions [12], (jjfi.o by using the eigenvalue equation (10) obtained 
in the previous section. Writing A = A^, + iXi and separating Eq. (10) into its real and imaginary parts we get 



Since the perturbations u{x,t) corresponding to n = again yield synchronous solutions, the linear stability of the 
synchronous state requires that all solutions of Eq. (10) have Xfj < for all non-zero integer values of n. The marginal 
stability curve in the parameter space of (k, r^) is defined by Xji — and in principle can be obtained by setting 
Xr = in Eq. (14), solving it for A/ and substituting it in Eq. (15). In practice it is not possible to carry out such a 
procedure analytically for the integral Eqs. (14, 15) and one needs to adopt a numerical approach, which is discussed 
in the next section. 

To systematically determine the eigenvalues of Eq. (10) in a given region of the complex plane we use multiple 
methods in a complementary fashion. First, Eq. (11) is solved for fl for a given set of values of the parameters k, lo 
and T„i. Following that, we need to determine the complex zeros of the function /(A) defined as 



which is equivalent to finding solutions A of (10). To do this we have primarily relied on the numerical technique 
developed by Delves and Lyness [14] based on the Cauchy's argument principle. By this principle the number of 
unstable roots m of /(A) is given by 



where the closed contour C encloses a domain in the right half of the complex A plane with the imaginary axis forming 
its left boundary. Once we get a finite number for m we further trace the location of the roots by plotting the zero 
value contour lines of the real and imaginary parts of the function /(A) in a finite region of the complex plane (A/?, A/). 
The intersections of the two sets of contours locate all the eigenvalues of Eq. (10) in the given region of the complex 
plane. The computations are done on a fine enough grid (typically 80 x 80) to get a good resolution. A systematic 
scan for unstable roots is made by repeating the above procedure for many values of the perturbation number n and 
by gradually extending the region of the complex plane. We have made extensive use of Mathematica in obtaining 
the numerical results on the stability of the synchronous states. 

In Fig. 1 the solid portion of the curve shows the stable synchronous states of Eq. (3) for k = 2.0 and for various 
values of w and f. The terminal point on a given solid curve of fl vs f marks the marginal stability point. The 
marginal stability point is seen to shift towards larger values of f as one moves down to curves with lower values of oj. 
A more compact representation is obtained if one plots rt — u versus ftf, since in this case the solutions corresponding 
to different values of lu for a given k consolidate onto a single curve, as shown in Fig. 2 for k = 0.05, 2.0 and 10.0 
respectively. 

It is seen from both figures that the stability domains of the synchronous solutions are restricted to the lowest 
branch where the curves are decreasing. This suggests a heuristic necessary (but not sufficient) condition for stability 
of the synchronous solutions: From Fig. 1 we have that dVl/df < for stable synchronous solutions, and from Eq. (11) 
we calculate 




(14) 



(15) 
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where 




(17) 



and 



K(e'" - 1) 
e" — 1 — K 
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fir 

FIG. 2. Solutions of Eq. (11) for the synchronous oscillation frequency for several values of k, plotted in terms of f2 — a; versus 
Qf. Note that Q, — uj = H{Q.f,Hi) for the equilibrium solutions (see Eq.(19)). In this representation, the different curves of 
Fig.l corresponding to different values of u collapse onto a single curve for a given value of k. The solid (black) portions of the 
curves correspond to stable synchronous states and the dotted ones (red) to unstable synchronous states. 

Since > and / is bounded, the denominator in Eq. (16) is positive for small values of f . Hence, for positive fJ, 
the requirement dVl/df < implies the condition J > 0, that is, 

J \z\G{z) cos {c^nf\z\)dz > (18) 

An alternative approach to arrive at the necessary condition (Eq.l8) would be to make use of the results presented 
in Fig. 2. We recast the dispersion relation given by Eq.(ll) in the form: 

n^uj = H{nf,K) (19) 

where 

H{Qf, k) — — J G{z) sin (c„ilf |2:|) dz 

It is seen from Fig. 2 that the stability domain of the synchronous solutions is restricted to the lowest branch where 
the curves have a negative slope. This again suggests a heuristic necessary condition for stability of the synchronous 
solutions to he H' < leading to the necessary condition given by Eq.(18). The prime indicates a derivative of 
H{flf,K,) w.r.t flf. 

As we will see later, the marginal stability curve obtained from H' or I — does lie above the true marginal stability 
curve (see Fig. 4), confirming that condition (18) is necessary but not sufficient for the stability of synchronous states. 

The points where solid and dotted lines meet in the curves of Fig. 2 mark the marginal stability point for the 
respective k values. These points are obtained for a range of k values and are plotted in (fir, k) space in Fig. 3 by 
filled points. They all lie on a single curve, which is analytically derived below. Our numerical results further reveal 
that for the marginal stability points, the imaginary part of the eigenvalue of the mode is zero — in other words the 
mode loses stability through a saddle-node bifurcation. It can easily be checked by inspection that A/ = is one 
of the solutions of Eq. (15) for any value of A^j; however, it is not evident analytically that this is the only possible 
solution for A^ = 0, and our numerical results have helped us confirm that this is indeed the case. Hence, putting 
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FIG. 3. The marginal stability curve (solid curve) in the (Or, k) space, obtained from the lowest branch solutions of Eq. (20) 
for n — 1. The filled circles correspond to numerical results from eigenvalue analysis of Eq. (10), and show a perfect fit to 
the analytical result. The dashed and dotted curves correspond to marginal stability curves obtained for n — 2 and n — 3 
perturbations, respectively. The symbols S and U denote stable and unstable regions in the parameter space. 



I=H'=0 




FIG. 4. A numerical fit to the marginal stability curve gives an approximate scaling law in the form of an offset exponential 
relation between fir and k. The marginal stability curve (solid, in black) of Fig. 3 has been replotted along with the fitted 
curve (dashed, in blue). The dotted curve (in red) is obtained from the condition H' or I = 0. 
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A_R = A/ = in Eq. (14) we get the following integral relation between the parameters fJ, r™ and k. 

1 

G{z) cos {nT,n\z\) [1 - COS (nnz)] dz = (20) 

-1 

Further, we have also observed that the most unstable perturbation is the one with the lowest mode number, namely 
n = 1. Therefore Eq. (20) with n = 1 defines the marginal stability curve, so the condition for synchronization takes 
the form 

1 

G{z) cos {nT^\z\) [1 - cos (ttz)] dz>0 (21) 

1 

The solid line in Fig. 3 is the analytical curve of marginal stability defined by setting the left side of Eq. (21) to zero, 
and it can be seen that the numerically calculated marginal values (represented by points) fit this curve perfectly. 
The figure also shows the stability curves obtained for the n = 2 and n = 3 perturbations (dashed and dotted lines, 
respectively) and these are seen to lie above the n = 1 marginal stability curve. We have carried out a numerical 
check for a whole range of higher n numbers and the results are consistent with the above findings. 

For a system with constant delay r (i.e. if Tm\z\ is replaced by r in Eq. (3)), the cos(r2T) term can be taken outside 
of the integral in (21) and the remaining integrand is nonnegative. Hence, the synchronization condition in this case 
becomes simply 

cos(f7r) > 0. (22) 

This agrees with the results obtained previously for constant-delay systems [5, 15]. Thus our result, as given by Eq. (20), 
generalizes the condition (22) to systems with space-dependent delays, and shows a nontrivial relation between the 
spatial connectivity and delays for the latter case. 

In order to gain some intuition into the complex interaction between connectivity and delays, we have obtained an 
approximate expression for the marginal stability curve by a numerical fitting procedure, yielding the relation 

r^f < 0.58 + 0.56e-°-3'*'" (23) 

for the stability of synchronous oscillations. Here, the left side involves the temporal scales of the dynamics (namely, 
it is the average time delay normalized by the oscillation period of the synchronized solution) while the right hand 
side involves the spatial scales of connectivity. In this view, the synchronization condition is a balance between 
the temporal and spatial scales. For high connectivity {k 0), the system can tolerate higher average delays in 
maintaining synchrony, and the largest allowable delays decrease roughly exponentially as the spatial connectivity is 
decreased. In the same figure we have also plotted with dotted curve the approximate condition (18), which is found 
to lie above the marginal stability curve in the entire range of k. The disparity between the two curves becomes 
particularly noticeable at large values of k. 

STABILITY CONDITIONS IN THE LIMITING CASES 

Since the marginal stability curve given by the condition Eq. (21) is an analytical one, one can obtain the limiting 
values of the phase shifts (fif ) in the global {k — > 0) as well as in the local limits (k -^• oo). They are given as: 

Qf < = 1.11072 (24) 
2v 2 



in the case of global coupling and 



nf < ^= 0.57735 (25) 
v3 



in the case of local coupling. 

Similarly one can obtain the limits on the frequency depression for the stability of the corresponding synchronous 
state in the two limiting cases (see Fig. 2). In the case of global coupling the condition becomes: 



n-uj> = -0.722819 (26) 



2^2 
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whereas in the case of local coupling the condition is : 

- a; -0.4433013 (27) 

4 

We note from Fig. 2 that the disparity between and w (|f2 — increases with delay induced phase shifts (fir ) 
and at some critical point the synchronous state becomes unstable. The value of the disparity at this critical juncture 
is larger for higher connectivity (global) and smaller for lesser connectivity (local) as we see from above. 



CONCLUSIONS AND DISCUSSION 



We have investigated the existence and stability of the synchronous solutions of a continuum of nonlocally coupled 
phase oscillators with distance-dependent time delays. Our model system is a generalization of the original Kuramoto 
model by the inclusion of naturally occurring propagation delays. The existence regions of the equilibrium synchronous 
and traveling wave solutions of this system are delineated. The equilibrium solutions of the lowest branch are seen to 
exhibit frequency suppression as a function of the mean time delay. We have carried out a linear stability analysis of the 
synchronous solutions and obtained a comprehensive marginal stability curve in the parameter domain of the system. 
Our numerical results show that the synchronous states become unstable via a saddle-node bifurcation process and 
the most unstable perturbation corresponds to an n = 1 (or kink type) spatial perturbation on the ring of oscillators. 
These findings allow us to define an analytic relation, given by Eq. (21), as a condition for synchronization. We have 
also obtained approximate forms for the synchronization condition that provides a convenient means of assessing the 
stability of synchronous states. Our results indicate an intricate relation between synchronization and connectivity 
in spatially extended systems with time delays. A detailed analysis of the stability of traveling wave solutions is 
currently in progress and will be reported later. 
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